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-  Summary 

•  i  ' '  .*•  '<  i-Ss;'} - 

-  We  discuss  a  method  for  the  estimation  of  unknown 
parameters  (variable  as  well  as  constant)  occurring  In 
a  hyperbolic  system.  In  the  context  of  a  seismic  appli¬ 
cation.  iip  present  both  theoretical  results  and  some 
numerical /(test)  examples.  ^ 


i  '  ‘  *  Introduction 

We  have  developed  a  numerical  algorithm,  and  a 
corresponding  convergence  theory  to  solve  a  one 
dimensional  "seismic"  Inverse  problem.  The  response 
in  certain  classes  of  seismic  experiments  can  be 
modeled  by  the  following  hyperbolic  partial  differen¬ 
tial  equation  with  associated  boundary  and  Initial 
conditions: 

2 

^1  (*)  -^-jf  “  )  t  >  0,  x  C  [0,1] 


Hf(t.O)  +  k^u(t.O)  -  s(t;k) 


J?  (t,l)  +  k2  fx  (t’l> 


u(O.x)  *  ♦(x) 

(0.x)  *  *(x). 

Here  x  represents  depth  below  the  surface  of  the  earth, 
u  represents  displacement,  q^(x)  Is  the  mass  density 

of  the  medium  (in  the  most  general  case  an  unknown), 
and  q2(x)  is  an  unknown  elastic  modulus.  The  boundary 

condition  at  the  surface  (x«0)  Is  an  elastic  boundary 
condition  involving  the  unknown  (negative  constant) 
k^,  and  an  unknown  source  term  s(t;k).  For  our  treat¬ 
ment  of  the  problem  it  is  not  necessary  to  assume  that 
)s  Is  an  Impulse.  In  the  numerical  examples  presented 
below.  It  has  been  assumed  that  s(t;K)  has  a  known  k 
form,  and  only  the  unknown  k  (a  constant  vector  In  IT) 
Is  to  be  Identified,  but  this  also  Is  not  essential. 

We  show  for  example,  that  q2(x)  (similar  remarks  are 

valid  for  s)  can  be  identified  as  a  function  without 
a  priori  knowledge  of  its  shape.  The  Ideas  in  this 
case  are  similar  to  those  In  [2]  where  problems  with 
coefficients  which  are  unknown  functions  of  both  space 
and  time  are  discussed  for  parabolic  equations.  At  the 
bottom  boundary  (x«1),  an  absorbing  boundary  condition 
Is  Imposed,  involving  an  unknown  (positive  constant) 


The  purpose  of  this  condition  Is  to  limit  the 


Interval  of  computation  without  producing  artificial 
reflections;  it  allows  down-going  waves  to  pass  through 
the  boundary  undisturbed,  while  annihilating  up-going 
waves . 

We  assume  we  have  data  observations,  y^,  corres¬ 
ponding  to  uft^.O),  a  solution  of  (1)  evaluated  at  the 

surface.  The  Inverse,  or  Identification  problem  con¬ 
sists  of  minimizing  a  least-squares  function 
n  -  2 

J(q)  =  7  |y,  -  u(t4,0;q)|  over  an  appropriately  chosen 
1-1  1  1 

constraint  set  Q.  Here  (t.x)  ♦  u(t,x;q)  Is  the 
solution  of  (1)  corresponding  to  q(x)  s  (q^(x),  ci.,(x) , 

kj,k2>k).  We  follow  the  general  approach  developed  In 

[3]  and  [1];  we  first  reformulate  the  Identification 
problem  In  an  abstract  setting,  then  define  a  sequence 
of  approximate  finite  dimensional  Identification  prob¬ 
lems,  the  solution  of  vrfilch  generate  parameter  esti¬ 
mates  which  converge  to  a  solution  of  the  original 
Identification  problem. 

Convergence 

Motivated  by  the  fact  that  our  differential  equa¬ 
tion  can  be  written  as  a  system  using  the  variables 

(u,ut),  we  define  a  Hilbert  space  X(q)  s  V(q)  *  l2(q) 

where  V(q)  is  H1  (0,1 )  with  Inner  product 

1  2 
<v.w>v(q)  s  /  R2DvDwdx  -  q2(0)kjv(0)w(0),  and  L  (q)  Is 

0  °  1 
H  (0,1)  with  Inner  product  <v,w>p  q  5  ^  q^vwdx.  The 

X(q)  Inner  product  Is  then  given  by  .y^vjq) 

♦  <x2,y2>o,q  where  x  •  (xrx2)T,  y  ■  (yry2)T.  After 

a  straightforward  transformation  to  a  system  with 
homogeneous  boundary  conditions,  system  (1)  can  be 
rewritten  In  X(q)  as 

z(t)  ■  A(q)z(t)  ♦  G(t;q) 


where  z(t)€  x(q)  Is  Identified  with 


/u(t..)  \ 
''ut(t,.)/  * 


84  08 


boundary  conditions  are  Incorporated  Into  the  domain 
of  the  operator  A(q)  by  defining  Vg(q)  * 

(v  €  V(q)  O  H2(0,1)  |  Dv(0)  ♦  k,v(0)  -  0}  and 

domA(q)  ■  J(y)  ®  Vg(q)  *  H1  (0.1 )  |  v(l )  ♦  k20u(l )  -  oj, 
and  A(q)  Is  the  unbounded  linear  operator  defined  by 

Approved  for  public  releaso  { 

CO  distribution  unlimited. 


050 


V^n/q,  )0(q20)  0  J 

With  X (q )  and  A(q )  so  chosen,  for  each  q  A(q)  Is  a  dis¬ 
sipative  operator  in  X (q )  and  it  can  be  shown  that  in 
fact  Afq)  is  the  generator  of  a  Cg-semigroup  T(t;q)  on 

X (q ) .  Standard  semigroup  theory  can  then  be  used  to 
show  that  equation  (2)  has  a  unique  mild  solution: 

t 

i(t;q)  -  T(t;q)zn(q)  +  /  T(t-s;q)G(s;q)ds  (3) 

u  0 

and  the  identification  problem  can  be  restated  as 

n  i.  ,2 

(ID)  Minimize  J(q)  s  X  (t1 )  lx»0j  over  q  e  Q 

subject  to  z(sq)  satisfying  (3),  where  z^  repre¬ 
sents  the  first  component  of  z. 

Before  formulating  the  approximate  identification 
problems,  we  first  define  finite  dimensional  subspaces 

XN(q).  let  S3(aN)  be  the  subspace  of  C2  cubic  splines 
(as  in  [6],  pp.  78-81)  corresponding  to  the  partition 
&N  *  (x^.g  .  x1  *  i/N,  and  then  define  XN(q)  to  be 
that  subspace  of  S3(aN)  x  S3(aN)  which  satisfies  the 
boundary  conditions  corresponding  to  q,  i.e. , 

XN(q)  C  domA(q).  The  space  XN(q)  can  be  expressed  as 
the  span  of  a  set  of  2N  +  3  basis  elements,  which  are 

straightforward  modifications  of  the  standard  spline 

basis  elements  of  S3UN)  *  S3(aN)  (see  [5]  or  [4,  p.  38] 
for  details).  As  a  result  of  these  modifications,  the 
new  basis  elements,  and  thus  the  subspaces,  depend  on 
the  unknown  parameter  q.  It  is  clear  then, that  as  we 
Iterate  on  q,  these  spaces  will  change. 

One  assumption  we  make  about  the  constraint  set  0 
is  that  each  component  is  uniformly  bounded  above  and 
below,  implying  that  as  q  ranges  over  Q,  the  X(q)  norms 
will  be  uniformly  equivalent,  and  hence  the  spaces  X(q) 

will  be  equal  as  sets.  With  this  in  mind,  let  P*(q)  : 

N  9 

X (q )  -  X  (q)  be  the  orthogonal  projection  of  X(q)  onto 

X  (q)  in  the  X(q)  norm  (for  a  precise  statement  of  this, 
one  should  introduce  the  canonical  isomorphism  which 

associates  elements  of  X(q)  with  those  in  the  equivalent 
space  X(q),  but  to  shorten  this  presentation,  we  will 

omit  such  notation);  whenever  q  and  q  are  the  same  the 
projection  will  be  written  as  PN(q ) .  Define  AN(q)  » 

P  (q)A(q)P  (q)  and  define  the  approximate  system  in 
XN(q)  as 


tN(t)  -  AN(q)zN(t)  ♦  PN(q)G(t;q) 


ZN(0)  -  PH(q)z0(q)  . 

The  operator  A*(q)  inherits  the  dissipativi ty  of  A(q), 
and  is  also  the  generator  of  a  Cg-semigroup,  T^(t;q)  on 

x(q).  Moreover,  we  can  establish  the  existence  and 
uniqueness  of  mild  solutions  to  (4)  and  write  them  as 


ZN(t;q)«TN(t;q)PN(q)zn(q)r/  TN(t-s;q)PN(q)G(s;q)ds. 

0  (SI 

We  now  pose  the  approximate  identification  problem  as 

ft  2 

(IDN)  Minimize  j"(q)  s  j  |;i-zj(t1)|x.0|  over 

q  €  Q  subject  to  zN(-;q)  satisfying  (5). 

It  is  Important  to  note  that  since  XN(q)  1$  a 
finite  dimensional  space,  (4)  is  in  fact  an  initial 
value  problem  for  a  system  of  ordinary  differential 
equations.  Furthermore,  due  to  the  nature  of  the 
spline  basis  functions,  this  system  possesses  desirable 
numerical  properties;  for  example,  the  matrix  repre- 

|J 

sentation  for  A  (q)  Is  sparse  and  banded.  We  will  dis¬ 
cuss  below  a  spline  representation  for  q2(x),  which 

makes  solving  (4)  even  more  tractable. 

||J 

The  identification  problems  (ID  )  can  be  solved 
using  1HSL  routines  (a  levenberg-Marquardt  algorithm 
for  the  optimization.  Gear's  method  to  solve  the  dif¬ 
ferential  equations)  with  a  sequence  of  parameter  esti- 

u 

mates  (q  }  thus  generated.  One  then  would  like  to 
verify  that  this  sequence,  or  some  subsequence  thereof, 
converges  to  a  solution,  q*,  of  the  original  problem, 
(10). 

We  use  the  following  version  of  the  Trotter-Kato 
Theorem. 

Theorem:  Let  (B,|*|)  and  (BN,|-|N),  N  «  1,2,...,  be 

Banach  spaces  and  let  *N  :  B  +  BN  be  bounded  linear 

Operators.  Further  assume  that  T(t)  and  TN(t)  are  C„- 
N  u 

semigroups  on  B  and  B  with  infinitesimal  generators 

A  and  AN,  respectively.  If 


(i)  11m  |„"x| 


for  all  x  £  B, 


(11)  there  exist  constants  M,  u  independent  of  N 
such  that  |TN(t)|N  £  Meut,  for  t  >  0, 

(iii)  there  exists  a  set  DC  B,  o  cdom(A),  with 

(Xq-A)d  «  B  for  some  xQ  >  0,  such  that  for  all 
x  £  D  we  have 

|aW  -  ttNAx  |  ^  -  0  as  N  -  » 

then  |TN(t)A  -  *NT(t)x|N  -  0  as  N  -  »,  for  all  x€  B, 
uniformly  in  t  on  compact  intervals  in  (0,»). 

We  apply  this  theorem  by  identifying  (B,|-|)  with 
(X(q),  |  •  | -),  (BN.MN)  **1th  (X(qN),|-|  N).  A  with  A(q), 
and  AN  with  AN(qN),  and  we  assume  (qN)t,is  an  arbitrary 
sequence  such  that  qM  *  q  in  an  appropriate  sense. 

Any  version  of  the  Trotter-Kato  Theorem  would 
Involve  some  convergence  statement  about  the  generators 
AN(qN)  and  A($).  As  mentioned  earlier,  the  spaces 
X  (q  )  and  the  domains  of  the  approximate  operators 
AN(qN)  change  with  qN.  Moreover,  elements  of  domAN(qN) 

U 

satisfy  the  boundary  conditions  corresponding  to  q  , 


b 
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while  elements  of  domA(q)_  satisfy  the  boundary  con¬ 
ditions  corresponding  to  q,  and  in  general  there  is  no 
inclusion  relation  between  these  sets.  This  necessi- 

u  ij 

tates  the  use  of  an  operator  *  :  X(q)  ♦  X(q  )  which 

maps  elements  of  domA(q)  into  those  of  domA(q  )  so 
that  it  will  be  possible  to  compare  these  elements. 

Once  the  Trotter-Kato  Theorem  has  been  used  to 
show  the  convergence  of  the  semigroups,  it  can  be  shown 

that  also  the  (mild)  solutions,  zN(t;qN)  of  (4)  con¬ 
verge  in  X(qN)  to  the  (mild)  solution  z(t;q)  of  (2) 
(again,  a  precise  statement  of  this  convergence  would 

require  the  use  of  the  canonical  isomorphism)  whenever 
N 

q  -  q  in  an  appropriate  sense.  With  this  result  and 
the  following  theorem  (from  [4]  or  [5])  it  can  be  shown 
that  q*  is  a  solution  to  the  inverse  problem. 

1  24k 

Theorem:  Assume  Q  is  compact  in  the  C  *  H 

topology.  If  q  -  zQ ( q ) ,  q  -  PN(q)z,  q  -  TN(t;q)z, 

z  £  X (q )  are  continuous  in  this  same  Q-topology,  with 
the  latter  uniformly  in  t  €  [0,T],  then 

N  N 

(i)  there  exists  for  each  N  a  solution  q  of  (ID  ) 

•  N 

and  the  sequence  (q  }  possesses  a  convergent 
subsequence  qN|t  ■*  q. 

(ii)  If  we  further  assume  that,  for  any  sequence 
{q3  1  in  Q  with  q3  ♦  q,  we  have  z3(t;q3)  -*■ 
z(t;q)  as  j  -  -  uniformly  in  t  €  £0,T],  then  q 
is  a  solution  of  (ID). 

The  proofs  and  details  of  all  the  results  stated  above 
can  be  found  in  [5]  and  are  variations  of  the  general 
framework  developed  in  [3]. 
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Estimation  of  Functional  Coefficients 


If  we  make  further  smoothness  assumptions  on  the 
variable  coefficients  to  be  estimated,  and  stronger 
compactness  assumptions  on  Q.  we  can  search  for  an 
approximation  to  each  of  these  coefficients  as  a 
finite  linear  combination  of  cubic  splines,  reducing 
the  infinite  dimensional  optimization  problem  to  a 
finite  dimensional  one.  In  our  numerical  examples,  we 
assume  for  computational  ease  that  q^(x)  =  1,  and 

search  for  q2(x)  in  a  function  space.  For  rotational 
convenience,  we  will  write  Q  as 
if  q 


Qj  *  Q3,  so  that 

(q^(x),q2(x),k.|,k2»  k )  £  Q,  then  q  ^  €  Q  ^ ,  ^2  ^  Q2  * 
and  (kpk^.^G  Q^.  If  we  let  IM( q2 ( * ) )  denote  the 
interpolate  of  q2  in  the  space  of  cubic  splines  S3(aM), 
then  following  arguments  similar  to  those  in  [2],  we 

can  conclude  that  I  (Q_)  has  the  following  represen¬ 
tation:  c 


q^*)  :  [0.1]- 


q^l*) 


M*3 


l  ciB1 (x),ci 


e  Ci 


i*i 


where  (B^(x))  are  the  basis  functions  for  S3(cM)  and 


each  c,  is  a  compact  subset  of  ft.  Searching  for  (an 

I  u 


KnH  approximate)  qj  in  IM(Q2)  1*  then  equivalent  to  search- 

<nn  fny  tr  r  i-  II,  r  .  f  .  r  r  ^ 


ing  for  (c^.Cj . c^j)  in  C1  *  C2 
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Define  QM  a  {(IN(q2).k1.k2.k)  |  q2  £  €  o]i. 

We  can  fix  M  and  iteratively  solve  (ID*)  over  QM  to 
obtain  a  sequence  {qN(M))N  with  q%)  s  (q^tMhkj.k*, 
k*)  and  q*(M)€  I*^).  The  convergence  theorems 

discussed  in  the  previous  section  guarantee  that  this 
sequence  contains  a  subsequence  (relabeled  for  con¬ 
venience)  such  that  q*(M)  ♦  q(M),  where  q(M)  satisfies 
J(q(M))  1  J(q)  Tor  any  q  €  Q That  this  sequence 
{q(M))M  is  In  fact  a  set  of  good  approximations  to  a 

solution  of  (ID)  can  be  argued  as  follows.  Under  the 

M 

proper  compactness  assumptions  about  Q,  Q  will  also  be 
compact,  ensuring  the  existence  of  a  subsequence 

(relabeled,  if  necessary  )  of  q(H)  such  that  q(M)  ♦  q 
in  Q.  Under  further  (although  not  too  restrictive) 
assumptions,  among  them  that  the  solution  z(t;q)  is 
continuous  in  q  (which  can  be  proved  using  the  Trotter- 
Kato  Theorem),  one  can  prove  that  q*  is  a  solution  to 
(ID).  This  proof  involves  the  compactness  of  Q  and 
standard  spline  error  estimates  such  as  those  found  in 

m. 


Numerical  Examples 

In  the  examples  to  be  presented  below,  the  "data" 
has  been  generated  using  an  Independent  finite  dif¬ 
ference  scheme,  urttere  known  "true"  values  of  the  param¬ 
eters  were  preassigned.  We  begin  each  example  with  an 

initial  guess  q®  and  a  choice  of  N  and  solve  (ID*)  to 

ij  u 

obtain  q  .  We  then  use  this  q  as  the  initial  guess 
for  the  next  value  of  N.  All  examples  were  produced 
either  on  an  IBM  VM/370  or  a  CDC  6600. 

Example  1:  For  this  example  we  parameterized  q2  as 
q2(x)  =  3/2  +  (l/»)  tan*1[q21(x-  q22)].  where  q21  and 
q22  are  to  be  estimated.  We  used  s(t;k)  ■  0.  and 

Initial  conditions  a(x)  *  ex,  ip (x )  *  -  3ex.  Data 
points  were  chosen  at  x  *  0  and  fifteen  equally  spaced 
time  values  in  (0,1].  We  obtained  the  following 
results: 


<21 

4 

kN 

*1 

kN 

*2 

uV) 

N*  4 

5.873 

0.503 

-0.995 

3.005 

0.15*  10'3 

N*  8 

5.929 

0.497 

-1.001 

3.001 

0.12*  10*4 

True 

Values 

6.0 

0.5 

-1.0 

3.0 

Start,, 
up  (q°)  ; 

5.0 

1.0 

( 

fo 

o 

1 

2.0 

Example  2:  We  added  random  noise  to  the  data  in  this 
example,  at  a  level  of  about  3*.  We  searched  for  q2 

as  a  constant,  we  used  s(t;k)  «  kj (l-e*5t)ek2t,  and 

used  zero  initial  conditions.  Data  points  were  chosen 
at  x  *  0  and  fifteen  equally  spaced  time  values  in 
(0,2].  We  obtained  these  results: 


N 

q2 

kN 

X1 

kN 

k2 

kH 

1 

kN 

*2 

jV) 

N*  4 

2.887 

2.076 

[  -0.996 

-2.110 

1.009 

0.12*10'3 

N  ■  8 

!  2.947 

2.032 

-1.013 

-2.038 

1.027 

0.17*10’3 

True 

Values 

3.0 

2.0 

-1.0 

1 

-2.0 

1.0 

Start- 
Up  (q°) 

2.0 

1.5 

!  -o  s 

-1.0 

2.0 

Exanplc  3:  In  this  example,  we  searched  for  q^(x)  In 
the  space  of  cubic  splines;  the  true  q2  used  was 
q2(x)  *  1.5  tanh[6(x -  •»)].  We  used  s(t;k)  *  0,  and 

♦(x)  *  ex,  e(x)  «  -  3ex.  We  did  not  search  for  the 
boundary  parameters  in  this  example;  the  true  values, 

*  * 

kj  ■  -1.0,  k2  •  3.0  were  used.  The  data  points  were 

chosen  at  seven  equally  spaced  spatial  values  In  [0,1] 
and  three  equally  spaced  time  values  In  (0,1].  Our 
Initial  guess  for  q2(x)  was  the  constant  function 

q2(x)  »  6.0.  With  N  *  4  (for  the  state  approximation) 
and  M  «  1  (coefficient  approximation)  we  obtained  an 
estimate,  q*.  for  our  functional  coefficient  such  that 
|q2  *  q2l0  *  0.099,  and  J4(q4)  «  0.48*  10'2.  We  have 

several  spatial  observations  in  this  example  rather 
than  only  the  one  at  the  surface;  this  Is  more  repre¬ 
sentative  of  problems  that  arise  In  treating  data  from 
"bore-hole"  type  of  seismic  experiments.  In  which 
receivers  are  located  at  various  points  down  a  well. 
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